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I. INTRODUCTION 

Strongly correlated electron systems have received a 
great deal of attention in the last twenty years, owing 
to the interest in classes of materials such as high tem- 
perature superconductors or heavy fcrmions. Some of 
the striking properties of these materials such as strong 
mass renormalizations, Mott insulating phases or uncon- 
ventional magnetic properties are clearly due to the cor- 
relation between electrons, an aspect ignored or poorly 
taken into account in conventional band theories. This 
has led to the development of an entirely new field and of 
new theoretical schemes and techniques. Among those, 
Dynamical Mean- Field Theory (DMFT)i has emerged 
as one of the most powerful, both for model Hamiltonians 
and as a way to take correlations into account in realistic 
electronic structure calculations 2 . 

Within DMFT, spatial correlations are frozen, while 
local quantum dynamics is fully preserved, as it happens 
in the infinite coordination limit, where DMFT becomes 
indeed the exact theory. Under this approximation, a lat- 
tice model finds an effective description in terms of an im- 
purity model in which an interacting site hybridizes with 
an effective bath of free electrons. The mapping onto the 
impurity model is enforced by a self-consistency condi- 
tion^ which contains the information about the original 
lattice. The self-consistency equation, as we will see, con- 
nects the hybridization function of the impurity model to 
the local Green's function. Therefore we can solve a lat- 
tice model within DMFT once we are provided with a 
method to solve the impurity model and compute the 
Green's function. The Anderson impurity model (AIM), 
albeit much easier to solve than the original lattice model, 
is still a non trivial many-body problem whose solution 
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requires either approximations or the use of numerical 
methods. Both "exact" numerical methods (Exact Diag- 
onalization (ED) 4 , Quantum Monte Carlo (QMC)£ and 
the recently introduced continuous-time version^, Nu- 
merical Renormalization Group (NRG )£,...) and approx- 
imate "analytical" methods (Iterated Perturbation The- 
ory^, the Non-Crossing Approximation^ and its slave- 
rotors extensions^, the self-energy functional method 10 , 
and others) have been successfully employed. However, 
most methods have limitations confining them to either 
a specific regime (e.g., high temperatures), or to the in- 
vestigation of specific physical aspects (e.g., low energy 
quantities). In particular, focusing on numerical meth- 
ods, Hirsch-Fye QMC is well suited for relatively high 
temperatures (and weak to intermediate correlations), 
while ED based on the Lanczos method has been up to 
now used only for T = 0. The NRG, which uses Wilson's 
scheme to solve the AIM is perfectly suited for an ex- 
tremely accurate determination of the low-energy part of 
the spectra at zero temperature, but it is slightly less ac- 
curate on the high-energy part of the spectrum, and for 
finite temperatures. There is no established reliable tool 
to deal with the regime of finite but very low tempera- 
ture, which is particularly relevant in correlated systems 
in which very small energy scales arise, leading to sub- 
tle effects (such as spectral weight transfers) when the 
temperature is turned on. The aim of this paper is to in- 
troduce a simple modification of the Lanczos strategy in 
order to treat the low-temperature regime accurately and 
with a reasonable numerical effort. We emphasize that 
our approach is different from the Finite- Temperature 
Lanczos method developed by Jaklic and Prelovse k 11 ' 12 , 
which is built as a tool to use ED at any temperature, 
but it is in principle exact only at T — and in the 
large temperature limit. Our method is instead designed 
to treat the very low-temperature regime with the same 
accuracy of T = 0, while it can not be pushed beyond 
some mo del- dependent temperature without spoiling the 
rapidity of the Lanczos algorithm and thus without al- 
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most recovering the computational heaviness of the full 
diagonalization of the Hamiltonian. 

As we anticipated in the introduction, the DMFT 
method maps a lattice model onto an effective impurity 
model that we can write as 

-HaIM = ^2 e la a la a l° + Vl°(fl a l° + a L/o-) +^at (1) 
la la 

In this expression, /J and a, are creation operators for 
fermions in with spin a associated with the impurity site 
and with the state I of the effective bath, respectively. 
For simplicity we consider a single band model, but the 
formalism is easily extended to multiorbital models. -Hat 
is the on-site (atomic) part of the original lattice Hamil- 
tonian, which contains the interaction terms. For the 
Hubbard model, iJ at = £/ (n^ + ) + Un^nf, and © 
is an Anderson impurity model (AIM). A fundamental 
quantity is the so called dynamical Weiss field Qq [w), 
which describes the non-interacting part of the effective 
model, and it is related to the Anderson parameters Vi a 
and £i a by the relation 

Qq 1 (iu n ) =iu> n + M - ^ ■ Vl — • (2) 

Introducing the impurity Green's function G(r) — 
-(T t c(t)J(0)}, and its imaginary-frequency Fourier 
transform G(iu> n ), we can extract the impurity self- 
energy 

S(iWn) = So* fan) - G^fan)) (3) 

which within DMFT coincides with the local component 
of the lattice self-energy. 

The self-consistency equation which establishes the 
equivalence between the lattice and the impurity mod- 
els depends on the noninteracting density of states D(e) 
of the original lattice 

Gfan) = I de . D{£) (4) 

J lOJ n + fJL — e- 2^(lU>n) 

For an infinite- coordination Bcthc lattice with semi- 
circular density of states of bandwidth 2D, @ reads: 

D 2 

Q (iu) n ) — iu n + fj. — — G(iu n ). (5) 

A practical solution of DMFT consists of an iterative 
solution of the impurity model. Starting from a given 
choice of the Weiss field, the impurity Green's function 
has to be computed with some "impurity solver" . The 
knowledge of G allows to compute £ from which, exploit- 
ing the self-consistency condition J3J one finds a new 
Weiss field. The process is then iterated until conver- 
gence. 

Let us now briefly recall the basic idea behind using 
ED as an impurity solver in the DMFT context. ED re- 
quires a truncation of the sums over / in Eqs. {1} and (JSJ) 



up to a finite value N s , the exact hybridization function 
being recovered in the limit N s — > oo. The accuracy of 
this method thus relies on how closely one can reproduce 
an infinite- N s bath with a finite- N s one. To be concrete, 
our discretized impurity model reads exactly as Eq. I|T|). 
with a small value of Ni. We can view this as a finite 
number of "sites" , each directly hybridized with the im- 
purity, in the so-called "star" geometry^ At every itera- 
tion, once Qq 1 is obtained through the self-consistence 
equation, the new set of Anderson's parameter is ob- 
tained through a fitting procedure, where a functional 
distance between the Qq 1 coming from Eq. ([5]) and a 
discretized version is minimized. In this work we mini- 
mize the function 

X = W fan)\Gofan) ~ Go'fan)\, (6) 

n 

where Q N " is the inverse of the discretized version of the 
Weiss field, the norm | . . . | is the square root of sum of 
the squares of the differences of the real and imaginary 
parts, and W(iuj n ) is a weight function. In this work we 
take the flat function W(iu> n ) — 1, but more selective 
functions can be useful for specific problems. For exam- 
ple, one can give more weight to small frequencies using 
W[iu> n ) — 1 jui n ^ The truncation error measured by x is 
the only systematic error in the ED solution of DMFT. As 
shown in Refj», x decreases exponentially by increasing 
the number of levels N s , so that relatively small numbers 
provide accurate information. The method is also able 
to provide real-frequency quantities without the need of 
analytical continuation tools, even if the spectra are nec- 
essarily discrete, due to the discreteness of the effective 
model. Nevertheless, many informations about single- 
particle and optical spectra can be obtained, mainly as 
far as the evolution of spectral weight is concerned^ 

In order to access finite temperature properties, one 
needs in principle the full spectrum of the system. There- 
fore, the size of the matrix to be diagonalized (4 JVs+1 x 
4 Ns+1 ) poses severe limitations on the values of N s which 
can be handled. Even using all the symmetries of the 
Hamiltonian, one can hardly go beyond N s — 5,6 us- 
ing full diagonalization. According to when self- 
consistency is achieved, the hybridization function of the 
bath is proportional to the local Green's function for 
a Bethe lattice. Therefore, a rough approximation of 
the bath for small N s is equivalent to a poor descrip- 
tion of G(iuj) 16 . This is expected to be more relevant 
at low temperature, where the Green's function is more 
structured, while at high temperature some structures 
can be broadened and eventually washed out. Notice 
that, nonetheless, the value of the energies of the An- 
derson impurity model can adjust, and in particular, 
their value can become arbitrarily small, as well as the 
weights can vanish. In this way, the method is able to de- 
scribe, e.g., the Mott transition, where the Fermi- liquid 
coherence energy goes to zero. It is therefore desirable 
to increase N s in order to obtain a reliable description 
of the low-temperature region. A standard way to in- 
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crease N s is to replace a full diagonalization with the 
Lanczos algorithmic. In this method one builds an or- 
thonormal basis in the subspace spanned by the vectors 
| <f>) , H 1 4>) , H 2 1 0) , . . . , H Nl | (f)) , where | (f>) is an arbitrary ini- 
tial state with non-zero overlap to the groundstate. One 
can see that in this basis the Hamiltonian becomes tridi- 
agonal and that even severely truncating the basis (i.e. 
limiting the number of states in the Lanczos basis Ni) the 
lowest lying states converge to the exact ones very quickly 
as a function of jVj. In practice the groundstate is very 
well converged already for basis of the order of Ni ~ 100 
even for huge matrices of size of the order of millions. 
Further increasing Ni , the low-lying excited states gradu- 
ally converge with a speed which basically depends on the 
energy distance between those states. More care has to 
be taken to properly handle degenerate states and their 
multiplicity, as we briefly discuss in Sec. II Bl 

Due to these convergence properties, this method has 
been up to now used mainly for the investigation of zero- 
temperature properties, for which only the groundstate 
vector needs to be determined. Nevertheless, Lanczos 
diagonalization can in principle still be used at finite (but 
low enough) temperatures, for which just a few low- lying 
states are needed to describe the system. In this work, 
we demonstrate that such an extension of the Lanczos 
algorithm can be used successfully in the DMFT context. 
The next subsection describes our approach. 



A. Extension to Finite Temperature 

As we anticipated previously we present here a rather 
straightforward extension of the Lanczos scheme to finite 
temperature, in which a small number of excited states 
of the Hamiltonian matrix are computed. 

We first show that the relevant quantities of the 
AIM can be expressed as a sum over the eigenstates of 
the Hamiltonian, each with a Boltzmann factor which 
weights it according to the energy distance from the 
groundstate. Thus the sums can be truncated to a finite 
number at low- enough temperatures. This observation 
is trivial for the partition function Z = ^ n e~^ E " . For 
the impurity Green's function, we start from the usual 
spectral representation: 



G a {iu) n ) = ^=Y^ 



-0E n 



(7) 



in which \n) and E n are eigenvectors and eigenvalues of 
Haim- This can be recast in the form: 



G CT (iw n ) = i^e-^GS m) (^„) 



(8) 



where: 



E m — E n — iu) n E n — E m — iu) n 

n 

(9) 



The "partial" Green's function Ga(ioj n ) involves cre- 
ating (or destroying) a particle into state |m). It can 
be readily calculated from m) alone, without knowing 
the whole spectrum spanned by \n), just like the T = 
Green's function (which is simply (i^n)) can be com- 
puted from the ground-state^. The exponential factor 
in Eq. ([5]) indicates that at large enough j3 only a small 
number of eigenstates needs to be calculated, and the 
sums can be limited up to n — Nk ep t ■ Every other spec- 
tral quantity, like, e.g., the dynamical spin susceptibility 
can be cast in an analogous form exploiting the Lehmann 
spectral representation. 

The calculation of a few excited states (and hence in- 
vestigating very low temperatures) is obviously possible, 
even if it is not as straightforward as the evaluation of 
the pure groundstate. What is far less obvious is whether 
a reasonably manageable number of states is enough to 
access the temperature range in which the physical prop- 
erties of the system start to deviate significantly from 
T = (e.g., reaching the Fermi liquid coherence scale 
in the correlated metal). The answer to this question 
depends on the model and on the range of parameters, 
since it is mainly connected with the level spacing in each 
subsector with given quantum numbers. In this work, we 
address this question using as a benchmark test the half- 
filled Hubbard model. This amounts to iteratively solve 
the AIM Jl} , computing the Green's function, and deter- 
mine from it a new set of parameters ei , Vi by minimizing 
the difference between the two members of Eq. ([5]) . Then 
the new AIM is solved and the procedure is iterated until 
convergence. 

We demonstrate that the finite-T Lanczos procedure 
can be applied quite successfully to the investigation of 
the Mott transition region at finite temperature. Our 
main results are that: (i) The method allows us to eas- 
ily solve the model for N s = 6 at a considerably lower 
computational cost than full diagonalization (ii) Within 
finite-T Lanczos, larger values of N s , (N s =8,9 and in 
principle the same values that are accessible at T = 0) 
can be used, which require a huge computational effort 
using the standard full ED, where the Hamiltonian is 
fully diagonalized. (Hi) Using N s = 8 we can draw the 
phase diagram of the Mott transition up to temperatures 
close to the Mott transition point at a reasonable com- 
putational cost (i.e., keeping a relatively small number of 
states). 

We finally briefly comment on the difference between 
our use at finite temperature of the standard Lanczos 
algorithm and the well established finite-temperature 
Lanczos method developed by Jaklic and Prelovsekii. 
The method discussed in Refs. [TUl is an ingenious 
modification of the Lanczos algorithm, in which the ther- 
mal averages are obtained as averages over random sam- 
ples of shortened Lanczos chains. The original version 
of the method 11 provides remarkably good results in the 
relatively high-temperature regime, but it is not partic- 
ularly efficient at low temperatures, and modifications 
have been proposed to overcome this limitation 12 . On 



the other hand, our method is precisely built to provide 
basically exact results for low temperature, and it cer- 
tainly breaks down (or becomes infeasible) at some tem- 
perature. Thus, the two approaches are basically com- 
plementary. 



B. Control of the Approximation 

Since the main limitation of the Lanczos method comes 
from memory requirements, our finitc-T implementation 
can in principle solve matrices of the same size as for 
T = 0. In practice the evaluation of excited states nat- 
urally slows down the method, ultimately limiting the 
number of states we can handle. Notice that the com- 
putation of excited states does not only require a larger 
number of Lanczos steps, but it is further plagued by a 
loss of orthogonality in the Lanczos basis, which gives 
rise to the so-called "ghost states", i.e., to replicas of 
the converged vectors with small weight. Different pro- 
cedures have been devised to overcome this problem, 
mainly based on selective reorthogonalization 17 . This 
necessary complication of the algorithm leads to an in- 
crease of computational time, which depends on many 
details of the spectrum. 

We finally mention a potential limitation of the present 
approach, which descends from the relative capability of 
the Lanczos method to handle degenerate states. It is not 
difficult to realize that, if the matrix we try to diagonalizc 
has a degenerate spectrum, the algorithm is not able to 
separate the different states, and to properly determine 
the multiplicity. The simplest way to avoid this problem 
is to implement all the symmetries of the Hamiltonian, 
and to diagonalize independently the Hamiltonian matrix 
in each symmetry subsector. In this case all the degen- 
eracies associated to those symmetries can not plague 
the calculation, as the degenerate states will appear in 
separated subsectors. Therefore only the errors arising 
from accidental degeneracies can affect the accuracy of 
our calculation. At least in the case we discuss here, we 
hardly encounter any measurable effect of such degenera- 
cies, even if we can not completely rule out such effects 
in other models. 

Our approach has two sources of error, whose effects 
can be minimized in a partially conflicting way. The first 
is the standard discretization of the Weiss fields, mea- 
sured by the value of the distance x defined in which 
is more relevant for low temperature, and the second is 
the truncation in Eq. |[7J), which will obviously be more 
and more relevant as T is increased. 

We start by discussing the effect of the second kind 
of truncation, since the first has been already discussed 
in the literature, and it has been shown to be rather 
benign-. We show results for the paramagnetic half-filled 
Hubbard model, and we compare G{iuj n ) (we drop hence- 
forth the spin index) obtained from full diagonalization 
of the Hamiltonian matrix for N s = 6, which is basically 
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FIG. 1: (Color online) Imaginary part of the local Green's 
function on the Matsubara axis for different values of A^ ept 
compared with the full diagonalization of the Hamiltonian 
matrix for N s = 6, j3 = 50/D. The left panel shows U = 2 AD 
and the right one U = 2D. 



the maximum size which can be treated with full ED, 
with our results for different values of Nkept, for a rela- 
tively high temperature T = 1/(3 = 1/50. It is important 
to underline that Nkept is the total number of states in 
the full Hilbert space. Obviously the states will belong to 
different subsectors with given quantum numbers. This 
means that in each subsector the number of converged 
states will be significantly smaller, hence the calculation 
relatively fast. The comparison, reported in Fig.[T]clearly 
shows that the convergence of our method as a function 
of Nkept is extremely fast, and the results are indistin- 
guishable from the exact ones already for Nf. e pt = 10—20. 
Lowering the number of kept states, the result approaches 
the T = one. We notice that the results converge fast 
to the exact ones irrespective of the value of the inter- 
action, even in a "difficult" case such as U = 2. AD, for 
which the T = solution is a metal while the system is 
insulating at T = 1/50 (convergence is much smoother 
at e.g U = 2.0D). The inclusion of a few excited states 
is therefore enough to qualitatively modify the physics of 
the system. It is important to emphasize the significantly 
lower computational time of our method, in comparison 
to the full diagonalization. In our implementation of the 
selective re-orthogonalization, we gain a factor of ~ 10 for 
Nkept — 20 and ~ 20 for Nkept = 10 with respect to full 
diagonalization (The precise numbers depend on many 
details of the spectrum). In practice, the method only 
introduces a factor of around 3 for Nkept = 20 in the com- 
putational time with respect to T = ED, so it still sub- 
stantially faster than QMC methods. This benchmark 
of our approach allows us also to determine a criterion 
for stopping the inclusion of excited state. We define a 
"difference" introduced by the inclusion of the n th state, 
as D n = Y, iuin \G n - G n ~ 1 \ {G n being the Green's func- 
tion obtained by including Nkept = n states, i.e. n terms 
in the sum in Eq. ([5])) and stop when this distance be- 
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comes smaller than a given tolerance, whose value can be 
extracted from the comparison with full ED for N s < 6 
and exported to larger N s values where the full ED is not 
feasible. This is a first indication that perfectly affordable 
calculations provide essentially exact results from T = 
up to finite temperatures of physical interest. In partic- 
ular, we used our method for N s = 8, where the size of 
the Hilbert space inhibits, or makes extremely heavy, the 
full diagonalization of the Hamiltonian. The results, re- 
ported in Fig. [5J are definitely satisfactory. The Green's 
function obtained with our method for U — 2D, (3 = 60 
and Nkept = 40 is basically indistinguishable from QMC 
solution for the same physical parameters. 



merical value of U C 2(T) line by computing the local spin 
susceptibility, a quantity which dramatically changes at 
the transition point from a Pauli-likc susceptibility in the 
metal to a large (oc 1/T) value associated to local mo- 
ments in the insulator. This is physically related to the 
increase of the effective mass when the metallic behav- 
ior is lost. The characterization of U C ±(T) requires more 
care. While, like U C 2(T), this line is associated to the 
disappearance of a metastable solution, there is no ob- 
vious quantity with a critical behavior when this line is 
approached. In practice, at each /3, we moved from large 
to small U with extremely small steps, until the insu- 
lating solution disappears. Fig. [3] presents our results 
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FIG. 2: (Color online) Imaginary part of the local Green's 
function on the Matsubara axis from finite-T ED (N a = 8, 
jV fcept = 40) and Hirsch-Fye QMC 19 for (3 = 60 and U = 2D 



II. RESULTS 

In order to prove that our algorithm works in a wide 
range of parameters, we now draw a phase diagram for 
N s = 8 and Nkept = 40. This relatively high number 
of excited state has been chosen according to the cri- 
terion discussed previously. More precisely we obtain 
-D40 < 10~ 8 for the highest temperature we consider 
T = Q.02W, and obviously even smaller values for the 
lower temperatures. We notice that a larger number of 
states has to be used here with respect to N s — 6, due to 
the larger Hilbert space. The scenario for the Mott tran- 
sition in the paramagnetic sector in the Hubbard model 
is now well established. Two distinct solutions exist, with 
metallic and insulating character. The former exists for 
U smaller than a temperature-dependent value U C 2(T), 
and the latter for U > U c i(T). At T = the transition 
is of second order and takes place at U = U C 2(0), while it 
becomes of first order at finite temperature. The coexis- 
tence region U c \ < U < U C 2 shrinks as the temperature is 
increased and closes at a critical temperature T c , where 
the first- order line ends in a critical point. From a prac- 
tical point of view it turns out easier to determine the nu- 



FIG. 3: (Color online) Phase diagram for the Mott tran- 
sition in the paramagnetic sector obtained through finite- 
temperature ED with N s — 8 and Nkept = 40, compared with 
previous estimates by Numerical Renormalization Group of 
Ref. 20 and Quantum Monte Carlo The ED curves 

are stopped at the highest temperature where the chosen num- 
ber of states determined a negligible truncation error. 

for this phase diagram, and a comparison with the Nu- 
merical Renormalization Group results of Ref.— and the 
QMC results o&- and^, which are used as references 
of popular methods used to study the Mott transition. 
We did not add the results of other approaches (self- 
energy functional, continuous time Monte Carlo, . . . ) in 
order to make the figure more readable, and we empha- 
size that the aim of this comparison is to prove the ability 
of our approach to study finite-temperature properties 
accurately, rather than a detailed comparison with dif- 
ferent approaches. The reported data clearly show that 
our method not only reproduces the Mott transition sce- 
nario at a qualitative level, but also provides results in 
extremely good quantitative agreement with established 
methods. In particular our method is extremely close to 
NRG at low temperatures, where this approach is basi- 
cally exact, and it is in very good agreement with QMC 
at higher temperature, where the latter method becomes 
accurate. Our method therefore accurately bridges be- 
tween the most popular well established impurity solvers, 
and allows to span a sizable region of the phase diagram 
with good accuracy with a single approach. We find 
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that Nkept — 40 produces a negligible truncation error 
up to (3 = 50, where our solution is still in extremely 
good agreement with previous results. Unfortunately, it 
is apparently difficult to get closer to the Mott endpoint, 
where the number of states needed to get a reasonable 
accuracy becomes larger and larger, due to the critical 
fluctuations which tend to diverge as the critical point is 
approached. In principle one gets a Mott endpoint also 
with 40 states, but the large truncation error suggests us 
not to plot the data around this point, where the method 
becomes less reliable, at least quantitatively. 

A confirmation of the ability of our method to ac- 
curately describe the low-temperature regime, we cal- 
culated the inverse lifetime of the quasiparticles 1/t — 
Zqplirriuj^alrnY^iiuj), where Z qp = (1 — <9i?e£(a>) / '<9o>) _1 
is the quasiparticle weight. It has been shown that the 
metallic phase of the Hubbard model studied in DMFT 
is a Fermi-liquid. According to Landau Fermi-liquid the- 
ory, 1/t has to be proportional to T 2 at low tempera- 
tures, with a coefficient which increases as we approach 
the Mott transition. Our method correctly reproduces 
this behavior with limited computational effort, as shown 
in Fig. 2J This result is not easily accessible to standard 
methods, and it shows precisely the main virtue of our 
approach, which works at its best in the low-temperature 
regime, where Fermi-liquid behavior, and possible viola- 
tions are directly and unambiguously detectable. 
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FIG. 4: (Color online) Inverse quasiparticle lifetime 1/t for 
N s — 7 as a function of square temperature for the displayed 
values of the correlation strengths. The linear behavior in T 2 
characteristic of the Fermi liquid is apparent, with a slope 
that increases with U/W 



III. CONCLUSIONS 

We have shown that a finite- temperature extension of 
the Lanczos algorithm can be successfully applied to the 



solution of the self-consistent impurity model appearing 
in DMFT. The inclusion of a few excited states allows 
for a reliable description of the physics up to tempera- 
tures of the order of the Mott critical endpoint with a 
relatively small increase of computational cost with re- 
spect to the well-established T — exact diagonalization. 
The method allows for a computationally cheap investi- 
gation of single-band models with extreme accuracy at 
very low temperatures, where the Fermi-liquid behavior 
and its possible violations can be investigated unambigu- 
ously. Furthermore, the present approach opens the way 
to the use of ED as a (small) finite temperature solver 
for more timely lines of researches, like cluster extensions 
of the DMFT 23,24 , or realistic calculations of properties 
of correlated materials 2 , at basically the same computa- 
tional cost of T = studies. As we already discussed the 
present algorithm has indeed the same memory require- 
ments as the T = standard approach 



We notice that, despite the scaling of ED methods with 
the number of orbitals (or, equivalently, sites in the clus- 
ter) is not favorable, ED has been successfully applied 
at T = to three-orbital models^ and to CDMFT for a 
2x2 plaquette 2 ^. These implementations require a total 
number of levels of the order of N s ~ 12. As we discussed 
in Sec. II Bl the present approach can be applied to the 
same size of matrices (i.e., to the same values of N s ) ac- 
cessible to the T = method, and the only limitation 
is given by the computational time, that grows in order 
to obtain accurate excited states. The actual increase of 
total time will depend on the temperature and on the 
size of the system, but our results for N s = 8, where the 
increase factor is around 3 for a range of temperatures 
that approaches the Mott transition endpoint, are quite 
promising. We believe anyway that the main use of this 
approach for multiorbital or cluster models can be to elu- 
cidate the really small temperature range, which is never 
easy to capture with other impurity solvers in the DMFT 
framework. This range is reasonably accessible with an 
affordable increase of computational time for the largest 
systems used in CDMFT and multiorbital DMFT. 
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